function [theta_ML,success,fval_ML]=run_NAG(N_constr,lb,ub,a,theta0,user)

	ncnln=int64(0);
	nclin=int64(N_constr);
	nn=int64(size(theta0,1));
	mode=0;
	nstate=1;
	cjac=zeros(max(1,ncnln),length(size(theta0)));
	clamda=zeros(length(lb),1);
	r=zeros(length(theta0));
	ifail=-1;
	istate = zeros(size(lb,1), 1, 'int64');

	[cwsav,lwsav,iwsav,rwsav,ifail] = e04wb('e04uc');
	[lwsav, iwsav, rwsav, inform] = e04ue('Print Level = 1', lwsav, iwsav,rwsav);
	[lwsav, iwsav, rwsav, inform] = e04ue('nolist', lwsav, iwsav,rwsav);
	[lwsav, iwsav, rwsav, inform] = e04ue('Major Print Level==5', lwsav, iwsav,rwsav);
	[lwsav, iwsav, rwsav, inform] = e04ue('Col', lwsav, iwsav,rwsav);
	[lwsav, iwsav, rwsav, inform] = e04ue('derivative level=3', lwsav, iwsav,rwsav);
	[lwsav, iwsav, rwsav, inform] = e04ue('Hessian=no', lwsav, iwsav,rwsav);
	[lwsav, iwsav, rwsav, inform] = e04ue('Major Iteration Limit==6000', lwsav, iwsav,rwsav); 
	[lwsav, iwsav, rwsav, inform] = e04ue('Step Limit==2', lwsav, iwsav,rwsav); 
	[lwsav, iwsav, rwsav, inform] = e04ue('Minor Print Level==0', lwsav, iwsav,rwsav);
	[lwsav, iwsav, rwsav, inform] = e04ue('Monitoring File==-1', lwsav, iwsav,rwsav);
	[lwsav, iwsav, rwsav, inform] = e04ue('Verify Level == -1', lwsav, iwsav,rwsav);
	[lwsav, iwsav, rwsav, inform] = e04ue('Stop Objective Check at Variable == 46', lwsav, iwsav,rwsav);
	[lwsav, iwsav, rwsav, inform] = e04ue('Stop Constraint Check at Variable == 3', lwsav, iwsav,rwsav);
	[lwsav, iwsav, rwsav, inform] = e04ue('Function Precision == 1e-8', lwsav, iwsav,rwsav);
	[lwsav, iwsav, rwsav, inform] = e04ue('Optimality Tolerance == 1e-4', lwsav, iwsav,rwsav);

	[~, ~, ~, ~, ~, fval, ~, ~, x, ~, ~, ~,~, ifail] = e04uc_global(a, lb, ub, '', 'obj_MLE', istate, cjac, clamda,r, theta0, lwsav, iwsav, rwsav,'n', nn, 'nclin', nclin, 'ncnln',ncnln,'user',user);

	theta_ML=NaN*theta0(:,1);
	success=0;
	fval_ML=NaN;

	if ifail==0
		theta_ML=x;
		success=1;
		fval_ML=fval;
	end

	%
